# Code Summary

## Structure

The code computes time dynamics in three steps: First, the potential is computed for given aberrations for both internal rotational states. Second, the eigenstates of this potential are computed using a discrete variable representation (DVR), together with all coupling terms for the two-body case. Third, time evolution is computed. The file structure is as follows:

 1. Compute the aberration potential (Mathematica)
    1. `Jones.wl` Code library 1
    2. `Polarization.wl` Code library 2
    3. `compute-tweezer-potentials.m` Master file computing the potential and storing it as a h5 file.
 2. Compute the Hamiltonian (julia)
    1. `nacs-tweezer-dvr.jl` Contains some fundamental constants and the single-body Hamiltonian from a given potential.
    2. `nacs-compute-interactions.jl` Contains function to load the potential from the h5 file and functions to compute the two-body Hamiltonian in a basis of single-body eigenstates.
 3. Compute time evolution (julia)
    1. `master-equation-diffeqs.jl` Defines differential equations for a given Hamiltonian, to be used with julia's `DifferentialEquations` libraries
    2. `onebody-drive-dephasing.jl` Compute the one-body drive time evolution (Figure 3(d)).
    3. `onebody-echo-dephasing.jl` Compute the one-body spin-echo time evolution (Figure 3(c)).
    4. `onebody-xy8-dephasing.jl` Compute the one-body XY8 time evolution (Figure 3(a)).
    5. `params-onebody.jl` Contains the parameter for the onebody evolution and calls other necessary files.
    6. `nacs-spin-echo-twobody.jl` Compute the two-body evolution.
    7. `params-twobody.jl` Contains the parameter for the onebody evolution and calls other necessary files.

## Computing the potential

Note: For the code to work, all 3 files must be in the same directory.
Requires setting output path `of`.

The first part of the code defines units, fundamental constants, and experimental parameters (until line 96). The next part defines the simulation parameters, i.e. the aberrations and the system size, needed for convergence (lines 108 to 112).
The following part defines functions following the references to compute propagate the electric field.
The following part describes how to compute the potential for both internal states for a given electric field.
The following part finds the point of peak intensity (i.e. potential minimum for the ro-vibrational groundstate), and the best numerical integration method.
Finally, the potential for both internal states is computed on a grid around the peak and stored as h5 files.

## Hamiltonian from potential
These files contains functions, which are supposed to be called as illustrated e.g. in `nacs-spin-echo-twobody.jl`:
 - Lines 73-85 compute the single-particle Hamiltonian(s). 
 - Lines 86-107 compute the low-energy spectrum and eigenstates. Signs are adjusted to be mostly aligned, which is benefitial for two-body simulations.
 - `energy_offset` computes the code-detuning needed to match the physical detuning `detuning_pulse` for a pulse applied to molecules in their 3D motional ground state.
 - `energy_submanifold_arr` computes blocks in the Hamiltonian, if far-off resonant terms will be ignored, controlled by `energy_resolution`.
 - Lines 131-149 illustrate how to compute the two-body Hamiltonian in the presence of a pulse. The library `nacs-compute-interactions.jl` contains additional functions for the same task if some of the parameters in this function are zero. When switching to those confirm the correct order of left and right, which is only guaranteed for functions directly called in the file `nacs-spin-echo-twobody.jl`.
 - Analogous code for a single molecule can be found in `onebody-drive-dephasing.jl` lines 87-151.

A general scaling factor `trap_depth` is generelly used to scale the trap from 4MHz to `trap_depth`MHz.


## One-body simulation
We discuss the one-body at the example of `onebody-drive-dephasing.jl`. After line 151, first the liouvillians are computed. Time evolution for one time step is then given by their matrix exponential. In a loop, then all time steps are evolved. The results are saved in the directory `save_dir`. The other two files must be called after the `-drive` one, since they rely on pre-defined functions and Hamiltonians.

## Two-body simulation
During the pulses, no dephasing is included, and unitary evolution is given by exp(-i H t).
Analogous to the one-body simulation, we first define the dephasing operators and the initial thermal density matrix (`Vector` of disconnected blocks). Time evolution is computed using `KrylovKit.exponentiate`, which can handle time evolution of general arrays. Results are stored in `save_dir`.